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Solar Convection Simulations 
using a B-spline method 

By T. Hartlep and N. N. Mansour 



1. Motivation and objectives 

Convection plays a key role in energy transport and global circulation in the outer layers 
of the Sun, in generation of solar magnetic fields and in many phenomena associated with 
solar activity and variability. Observations of the solar surface reveal structures that have 
been classified as granules, mesogranules, and supergranules. The nomenclature reflects 
organization at three spatial scales ranging from about 1 Mm t o 30 Mm. 

Numerical simulations of the near surface region of the Sun ( Stein fe Nordlund 2000l) 
capture structures on the granular scale, but have not been able to detect organization 
at large scales. The physical mechanism of supergranulation is presently unknown. It 
has been suggested that supergranulation corresponds to a large convective cells which 
develop due to enhanced convective instability in the Hell ionization layer. This layer 
lies deeper than Stein & Nordlund have investigated so far. 

We know that as we go deeper from the surface of the sun, the turbulence structures 
become large, and the Mach number decreases. It is then advantageous to be able to 
change the spatial resolution in all three coordinates as a function of depth. In addition, it 
is numerically advantageous to use the an elastic approximation ( Qgura fc Phillips! 1962 ; 
Goughl [l969l : iGilman fc Glatzmaier 1981) to the compressible Navier-Stokes equations 



for deep domains. Using B-splines i n one coordinate d i rectio n and Fourier methods in 
the other two coordinate directions, Kravchenko et all (|l996l ) and Loulou et al. ( 1997 ) 
have developed B-spline- spectral numerical schem es that enable changing the resolution 
as a function of height. Kravchenko et al. ( 19961 ) applied the scheme to simulate the 
fully developed turbulent cha nnel flow by r e solvin g the near- wall region and relaxing the 
resolution in the core region. iLoulou et all (|l997[ ) simulated a fully developed turbulent 
pipe flow. They applied the scheme to remove the centerline singularity in cylindrical 
coordinates. 



2. Numerical method 

Under this effort, a new code to simulate a rectilinear section of the solar convection 
zone is being implement ed. As a first step, the numerica l method is implemented for a 



simple Boussinesq fluid (jOberbeck 118791 : iBoussinesd 1903) and is presented in this paper 
for this case only. Later, once the new code is fully functional and tested, equations based 
on the anelastic approximation will be used. 

2.1. Basic equations 

Using the temperature difference across the fluid layer, AT, the layer thickness d and the 
thermal diffusion time d 2 /k as units of temperature, length and time, the dimensionlcss 
Boussinesq equations read: 



d t v+ (v- V)u = -Vtt + PrV 2 v '+ RaPr0e z 



(2.1) 
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d t e+{v-V)9 = V 2 6 + v-e z , (2.2) 

V--u = 0. (2.3) 
v, 9 and ir denote the velocity of the fluid and the deviations of the temperature and 
pressure from their static profiles, i.e. the temperature and pressure profiles that would be 
observed in the case of no motion. The equations contain two dimensionless parameters: 
the Prandtl number Pr and the Rayleigh number Ra, which are defined as 

v „ gaATd 3 
Pr=~, Ra=- , (2.4) 

with k, ^, a and g being the coefficients of thermal conductivity, kinematic viscosity 
and thermal expansion, and the acceleration due to gravity, which points in negative 
z direction. The Boussinesq approximation assumes all these properties to be constant 
throughout the layer. 

2.2. Velocity decomposition 

A poloidal/toroidal decomposition is used for the velocity to automatically fulfill the 
continuity equation V • v = 0. Even though the divergence of the velocity does not vanish 
in the anelatic approximation, the equations that we are ultimately going to use, such a 
representation can still be used. Just not for the velocity itself, but for the quantity pov, 
where po denotes the reference density. 

In the Boussinesq case, the representation is in the form of 

v(x,y,z,t) = V x [ip(x,y,z,t)e z ] + V x V x [<j>{x, y, z, t)e z ] + U(z,t). (2.5) 

In this way the poloidal and toroidal scalars <p and ip are periodic in x and y, and U 
represents the horizontal average of the velocity, i.e 

U(z,t) = (v(x,y,z,t)) xy . (2.6) 

The z component of U must be zero to satisfy the continuity equation. U is thus a 
toroidal flow. However, the corresponding toroidal scalar varies linearly in x and y and 
is unbounded. Obviously, we require <f> to be periodic and to remain finite in the nu- 
merical representation. Th erefore U needs to be included in the decomposition of v 
( Schmitt k von Wah]|[l992h . 



Equations of motion for <j>, %p and U can be derived by evaluating the z component 
of the curl of the curl of (|2.ip , the z component of the curl of (|2.I|) and the horizontal 
average of (|2.I[) , respectively. 

2.3. Spectral method for the horizontal directions 

Periodic boundary conditions are used in the horizontal directions to reduce the influence 
of the horizontal boundaries on the flow. A natural choice is the use of Fourier modes in 
those directions, i.e. 



f(x, y, z,t) = J2 h& t)e-^ x+k ^ (/ = 0, V, 9), (2.7) 

k 



with k = k x e x + k y e y being the horizontal wave vector. Multiplying the equations of 
motion with weight functions e l ^ kxX+k y v ' > and integrating over the entire horizontal plane 
then yields equations for the z and t dependent Fourier coefficients: 



\n 2 



d t - Pr[dl - k 2 



4>i(^t)^n h {z,t), (2.8) 
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Figure 1 . Second-order B-splines denned on sets of 5 knot point in the case of equally spaced 
knot points (top diagram) and non-equally spaced knot points (bottom diagram). In both 
cases solid, dotted, dashed, dash-dotted, dash-dot-dotted and long-dashed lines denote splines 
Bi, . . . , B%, respectively. 



d t ~ Pr[d 2 z - k 2 } ^(z,t) = Kr(z,t), 



dt ~ [dl 



h(z,t)=K s (z,t). 



(2.9) 
(2.10) 



All non-linear terms are contained in the right-hand sides 1Z. The equation for U is very 
similar: 

[d t -Prd 2 }U(z,t) =Tl {z,t). (2.11) 

If fact, we can discuss the equation for ip^ 6? and U together, since they all can be 
written in the form 



[d t + C -C 2 d 2 }f(z,t) = K f (z,t) 
with some parameters Cq and C 2 - 



(2.12) 



2.4. B-spline method for the vertical direction 



Spatial discretization in the vertical direction z is done by a B-spline method (jKravchenko fc Moin 
Il998l) . i.e. the unknowns, 4>z,'iptiQz an d U, are expanded in terms of m-order piecewise 
polynomials called basis splines, or B-splines. With a given set of knot points {zrj, • • ■ , %n} 
that divide the z domain into N subintervals, one can construct N + m of these spline 
functions using the recursive expression 



B"\z) = 



— m — 1 



Zj—X Zj—wi—i 



-B r , 



A) 



1, 



,N + m). 

(2.13) 
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Here, B™ denotes the j'-th B-spline of order m. The 0-order splines are given by 

i o otherwise ' (2-14) 

As it turns out, m additional points at each side of the domain (z— m , . . . , z~i and 
zjv+ii • • • i ZN+m) are needed in the construction of B-splines near the boundaries. These 
virtual points are chosen to be equal to the boundary points zq and Zff, i.e. Z— m = . . . = 
z-i = zq and zjv+m = • • • = zjv+i = zjv- 

Second-order B-spline functions for two different sets of knot points are shown in 
figure Q] to illustrate how these functions look like. As can be seen from the figure, B- 
splines have compact support and are non-negative: 

r>m / \ J > if Zj- m -i < Z < Zj , , 

{Z) \ =0 otherwise ' [2Ab) 

and can be tuned to the specific needs at hand by the choice of knot points. In areas 
where small structures need to be resolve, e.g. in boundary layers, one can choose closely 
spaced knot points, while in areas where such high resolution is not required, the grid 
can be coarse to save computational costs. 

Incorporating a B-spline expansion of the unknowns, i.e. 

N+m 

f{z,t)= a fJ (t)B?(z) (/ = ^,^,^,l7), (2.16) 

into the governing equations ()2.8j) and (j2. 12[) . multiplying with weight functions Bf 1 and 
integrating over the z domain yields equations for the expansion coefficients atjit): 




(2.18) 

Three matrices, A4o, M.2 and A4^, occur in these equations, which contain integrals of 
the product of two B-splines and the product of a B-spline with a derivative of a B-splinc: 

(M n ) = B?8£B?dz (n = 0,2,4). (2.19) 

V / i=l,...,N+m J 

It follows from the property (12.15[) that all of these matrices are of band-diagonal form. 
Specifically, they have only m subdiagonals and m superdiagonals which contain non-zero 
entries. 

2.5. Time advancement 

A mixed implicit /explicit method using a Crank-Nicolson scheme for the diffusive terms 
and a second-order Adams-Bashforth scheme for the advection and buoyancy terms is 
used for the time advancement. 

In the end, one has to solve a matrix equation for each Fourier mode of each field 
(poloidal, toroidal, mean flow and temperature) in the form 



(2.20) 
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with df representing the vector of B-splinc coefficients at the new time step t n +i, i.e. 
Sf = (a f s(t n+1 ), . . . ,aj <N+m (t n+1 )) (/ again stands for either ^ ,0g or U), and 
the matrix A// is composed of the previously introduced matrices Mo, Mi and M^. TV/ 
is therefore band-diagonal and a fast method for solving band-diagonal systems can be 
used. 

2.6. Imposing boundary conditions 

So far, no boundary conditions have been imposed in the vertical direction. Since only 
one B-spline, and the derivative of only two B-splines are non-zero at the boundary, it 
is rather easy to implement various boundary conditions. If for example one needs to 
prescribe the value of f(z, t) at the boundary zo, e.g. f(z = zq, t) = g(t), one simply has 
to set the first B-spline coefficient to that value, i.e. a/ = g(t). The original equation 
for ctf^i has to be removed from the linear system (|2.20p , and at i's contribution has to 
be deducted from the right hand, i.e. equation (|2.20[) is modified in the following way: 
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More complicated boundary conditions in which the value of the B-spline coefficients at 
the boundary depend on the coefficients in the interior can be implemented by directly 
manipulating the matrix elements of Mf . 

2.7. Varying numerical resolution with depth 

The main advantage of such a B-spline method is the ability to vary the spatial resolution 
in all directions as a function of depth. As already mentioned, the choice of knot points 
in the construction of the B-spline functions determines the vertical resolution. Near the 
surface of the sun, where high spatial resolution is required, closely spaced knot points 
will be chosen, while the spacing can be increased in deeper parts of the convection zone. 
The spatial resolution in horizontal directions is determined by the number of Fourier 
modes that one wants to consider in the discretization: 

f(x, y,z,t)=J2J2 a f At)e- iik ' x+k " y) BT(z) (/ = ^, 6 S ), (2.21) 

o k 

and this number can be changed from one spline index j to the next. Since according 
to (|2.15| B-splines of different j overlap, the change in resolution occurs over several knot 
points. A consequence of the variation of horizontal resolution is, that for a given set of 
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wavenumbers k x and k y only a subset of B-splines has to be considered, namely those for 
which the maximum wavenumbers that we have chosen are larger or equal to k x and k y . 
The other coefficients are zero and the corresponding rows and columns in the matrix 
equations (|2.20p can therefor be dropped to reduce the numerical cost. 



3. Future work 

The implementation of the presented numerical method is nearly complete and will be 
tested shortly for some test cases of tur bulent Rayleigh -Benrad convection. Comparisons 
will be done with the results obtained bv lHartled l|2004h using a Chebyshev method. After 
that, the anelastic equations will be implemented in place of the Boussinesq equations. 

The final ingredient for the solar simulations is the definition of reasonable b oundary 
condit ions. As a start, we can use rather simple conditions like the ones used bv iMiesch 
(1998) and others: stress-free, impenetrable upper und lower boundaries with constant 
heat flux at the bottom and constant entropy at the top of the computational domain. A 
refinement to these boundary conditions would be the treatment of the upper surface as 
a free surface, which could be realized by additionally tracking a height function h(x, y, t) 
above the top end of the discretization (|2.21|) . 
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